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Abstract 

We present details of a new numerical code designed to study the formation and evapo- 
ration of 2-dimensional black holes within the CGHS model. We explain several elements of 
the scheme that are crucial to resolve the late-time behavior of the spacetime, including regu- 
larization of the field variables, compactification of the coordinates, the algebraic form of the 
discretized equations of motion, and the use of a modified Richardson extrapolation scheme to 
achieve high-order convergence. Physical interpretation of our results will be discussed in detail 
elsewhere. 
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1. Introduction 

The Callan-Giddings-Harvey-Strominger (CGHS) model [1] is a two-dimensional 
model of quantum gravity which has attracted attention due to the fact that it has black 
hole solutions with many of the qualitative features of four-dimensional black holes, 
while being technically easier to investigate. Various properties of black holes in this 
model, and other models inspired by it, have been studied extensively using analytical 
and numerical methods [2, 3, 4]; for pedagogical reviews see [5]. A recent focus point 
has been on using the CGHS model to investigate the black hole information loss prob- 
lem [6], where the importance of understanding the asymptotic behavior of the fields 
near right-future null infinity X\ was emphasized. In particular, sufficient conditions 
for the unitarity of the S-matrix were given. Although the full quantum equations are 
too complicated to solve, in the mean field approximation (MFA) the model reduces to 
a coupled set of non-linear partial differential equations, possessing a well-posed char- 
acteristic initial value formulation. Unfortunately, even for these equations, analytical 
solutions are not known except in special limiting cases. Therefore, to explore black 
hole formation and evaporation, numerical methods are essential. 

In this paper, we give details of the methods we have devised for accurate numer- 
ical calculations of the fields and related physical quantities in the CGHS model. We 
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give special attention to the macroscopic mass limit, which is the most challenging 
case to calculate and which has not been properly investigated before. An outline of 
the rest of the paper is as follows. In Sec. 2 we introduce the CGHS model, describe the 
variable definitions and conventions we use (which closely follows [6]), the analytical 
equations that we discretize, and the initial data we use. In Sec. 3, we describe some 
of the issues that would cause naive discretization of the equations to fail to uncover 
the full spacetime, and how to overcome them; this includes regularization of other- 
wise asymptotically-divergent field variables, compactification of the coordinates, the 
particular discretization scheme, and use of Richardson extrapolation ideas to increase 
the accuracy of the solution. In Sec. 3 we also discuss setting initial conditions near 
X, and how we extract the desired asymptotic properties of the solution. In Sec. 4 we 
describe various tests to demonstrate we have a stable, convergent numerical scheme 
to solve the CGHS equations. We summarize and conclude in Sec. 5. 

2. CGHS Model 

The action of the 2-dimensional (2D) CGHS model is given by 

S(g,4>,f) = y d 2 V e~ 2 * (R + Ag ab V a 4>V b 4> + 4« 2 ) 

- \j& 2 Vg ab V a fV b f . (1) 

where g ab is the metric, R is the Ricci scalar, is a dilaton field, / is a massless scalar 
field, G is Newton's constant and k is a constant of dimension inverse-length. Note 
that this action is similar to though not exactly the same as what would be obtained by 
dimensional reduction of the 4D Einstein-Klein-Gordon equations in spherical sym- 
metry. 

We are interested in metrics g ab that approach a Minkowski metric rj ab at past null 
infinity. We will denote the null coordinates of i] as z ± (see Fig. 1). Defining the fields 
$ and via 

$ = e- 2< ^ and g ab = O" 1 * ry ab = ft r/ ab , (2) 
we can write the equations of motion in terms of a set of evolution equations 

□ (9) / = □(,)/ = 

$<9+d_lne = -GT+_, (3) 

and constraint equations 

-d\ $ + d+ $<9+ In 9 = GT++ 

-d 2 _ $ + d- $<9_ In 6 = GT__, (4) 

where Or g ) (Dfr?)) is th e wave operator with respect to the metric g ab (r) a b), T ab is the 
scalar field stress-energy tensor with components denoted by T ++ , and we use the 
notation <9_ = d/d z -, and similarly for d + . Classically (and at the tree-level) T ab is 
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Figure 1: The Penrose diagram of the CGHS space-time at the mean field approximation (MFA) level for an 
incoming <5— function matter wave on null coordinates, showing the main features of the model. Unlike the 
3 + 1 dimensional case, here there are two past null infinities (left I£ and right X^ ) and two future null 
infinities (left X^ and right Xi), since the uncompact spatial coordinate is in the range (-co, oo), compared 
to a radial coordinate with range [0, oo) in 3 + 1 dimensions. The incoming <5— function matter wave forms 
a black hole. To the left of the matter wave, the space-time is exactly flat, and in the vicinity of X^ and X^ 
the space-time is asymptotically flat. The singularity and the dynamical horizon (dashed line) meet at finite 
z ± . The last ray is the null line connecting this point to X^. Within the MFA only the region of spacetime 
to the causal past of the last ray and singularity can (uniquely) be determined. 
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trace-free, hence T_| vanishes. However at the one-loop level T_| picks up a non- 
zero value due to the trace anomaly, which if we now consider the superposition of TV 
identical massless scalar fields is 

GT +- = ^ d + d - (InS-mO) (5) 

In a characteristic initial value problem, we specify initial data on a pair of inter- 
secting, null hypersurfaces z + (z~) = Zq and z~{z + ) = z^ , to the causal future of 
their intersection point (zq,Zq) (see [7] for a review). Thus one can see where the 
constraint equations (4) receive their name: for example, if we specify the scalar field 
/ (hence T ++ , T__) and metric field 9 on these surfaces as initial data, we are not 
free to choose which is then given by integrating (4). The constraint equations are 
propagated by the evolution equations (3), namely, if the constraints are satisfied on 
the initial hypersurfaces, solving for the fields to the causal future using (3) guarantees 
the constraints are satisfied for all time. This is exactly true at the analytical level, 
though in a numerical evolution this property of the field equations will in general only 
be satisfied to within the truncation error of the discretization scheme. 

Let {zq ,Zq) = {— oo, — oo ) be the initial data surface, and denote the choice of 
initial ingoing and outgoing scalar field profiles by 

/(*+*" = *o) = /+(*+) (6) 
f(z+ = z+,z-) = /_(*"). (7) 

where J + (zq) = /-(zq). Note that due to the conformally invariant nature of the 
wave operator in 2D (3), the solution for / over the entire spacetime is simply f(z + , z~) 
f+{z + ) + f-(z~) — /(zq, Zq). With these initial conditions, and the condition that 
the metric approaches Minkowski on left- and right-past null infinity, the solution to 
the constraints are [6]: 



e(^ ± ) = -k 2 x+x 

G 
2 



= 9(^)-^ dx+/; dx+(df+/dx+)i 



- jjfdx-ffdx-(df-/dx-)\ (8) 

where the notation F{z ± ) denotes evaluation of the corresponding field F on the given 
initial hypersurface, and 

K x+ = e KZ+ , KX - = -e- KZ ~. (9) 
For a first study, we will exclusively consider the case 

^(df + /dx+) 2 = M5(x+-x+), (10) 

with Xq = 1 (zq = 0), and no incoming matter from the left (/_ = 0). This choice 
reduces the problem to evolving the fields $ and according to (3) with the asymptotic 
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initial conditions 

e(^ ± ) = e K ( z+ -o 

$(z ± ) = e K{z+ - z ~ ] -M(e KZ+ - 1) , (11) 

for z + > 0, z~ -> -oo. Both fields are trivially given by e K ^ + ~ z ) for z + < 0. With 
these restrictions, any space-time is defined by the two quantities M and N. M is also 
the Bondi mass of this system as z~ — > — oo. 

The classical solution (h = 0) to the future of the delta-pulse matter wave is exactly 
given by (1 1), though now valid everywhere within this domain and not just the initial 
data surface. This spacetime contains an event horizon relative to right-null infinity, 
and a singularity inside (to the left of) the horizon. In the mean field approximation, 
the black hole evaporates, and the event horizon is replaced by a dynamical apparent 
horizon. It is expected that the full quantum theory will resolve the singularity, however 
in the MFA there is still a singularity inside the dynamical horizon. When the evapora- 
tion has proceeded to the point where the dynamical horizon meets the singularity (see 
Fig. 1), it becomes naked, i.e. visible to observers at Zj. The MFA equations cannot be 
solved beyond this Cauchy horizon, which we call the last ray. It should be possible to 
mathematically extend the spacetime beyond the last ray, in particular as the geometry 
does not appear to be singular here (except at the point the dynamical horizon meets 
the last ray), though we will not explore those issues here. 

In all our simulations we use G = h = k = 1. The CGHS model gives the same 
physics if N and M are scaled by the same number, and N/24 gives a natural scale 
for the unit mass. For the results presented here, we always use N = 24. Hence, by 
macroscopic mass, we mean M 1, and by sub-Planck-scale mass, we mean M<1. 

3. The Numerical Calculation 

In this section we describe several novel aspects of our solution scheme that al- 
lows us to uncover the physics of 2D black hole evaporation within the CGHS model. 
This includes compactification of the coordinates (Sec. 3.1), regularization of the fields 
(Sec. 3.2), the discretization and solution strategy (Sec. 3.3), and a Richardson extrap- 
olation algorithm to increase the order of convergence of the base, 2nd order accurate 
scheme (Sec. 3.4). We also discuss in Sec. 3.5 some difficulties in naively attempting 
to solve the discrete equations near null-infinity, and how we extract desired properties 
of the solution nearZ^ in Sec. 3.6. 

3.1. Compactification of the Coordinates 

Rather than discretizing the equations with respect to the z + , z~ coordinates, we 
introduce a compactified coordinate system z+ 6 [0, j] an d z c *= [0> !]• Use of com- 
pact coordinates is important for a couple of reasons, and essential for the M 3> 1 
case. First, to understand the asymptotic structure of the spacetime approaching Zj, 
it is useful to have the computational domain include Z^. Second, the uncompactified 
coordinate z~ is adapted to the flat metric near Z^; however, it turns out that most of 
the interesting features of black hole evaporation near the dynamical horizon occur in 
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Figure 2: A schematic view of the positions of the grid lines on the uncompactified space. Lines are concen- 
trated near the last ray, where we need higher resolution. They become distant as one approaches the null 
infinities. 

an exponentially small region Az~ ~ n~ 1 e~ GM > K before the last ray. One can think 
of this as essentially due to gravitational redshift. Classically (without evaporation), 
the redshift causes arbitrarily small lengths scales near the horizon to be expanded to 
large scales near 1^. Naively one might have expected that evaporation changes this 
pictures completely (as suggested by the Penrose diagram in Fig. 1). Instead, what we 
find is that although there is not an arbitrarily large redshift once back-reaction is in- 
cluded, there is still an exponential growth of scales, with the growth rate proportional 
to the mass of the black hole as indicated above. 

Thus, a uniform discretization in z~ that is able to resolve both the early dynamics 
near 1^, yet can adequately uncover the exponentially small scales (as measured in 
z~) of the late-time evaporation, will (for large M) result in a mesh too large to be able 
to solve the equations on using contemporary computer systems. To overcome this 
problem, we introduce a non-uniform compactification in z~ , schematically illustrated 
in Fig. 2), that provides sufficient resolution to resolve the spacetime near the last ray, 
yet does not over-resolve the region approaching . Specifically, the transformation 
from z~ to z~ we use is as follows. First, we relate the uncompactified z~ to an 
auxiliary (non-compact) coordinate z~ by 



where z~ G (— oo, z~ est ] and z~ G (— oo, 0]. z~ est is an estimate of the z~ coordinate 
of the last ray. This is also the earliest time in z~ that we will encounter the spacetime 
singularity, and at present we do not continue the computation past this point (the 
compactification functions can readily be adjusted to cover z~ G (— oo, oo) ). This 




(12) 
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way, the region near the last ray {z~ « z~ est , z~ « 0) is resolved by a factor of 
more than the regions away from the last ray. Next, we convert the auxiliary z~ to a 
compact coordinate z~ 

7- = _ e -Stan( 7r2c -- 7 r/2) + L j^- _ ^ _ (B) 

where S and L c are constants. This way, the last ray is located near z~ = 1. The 
relation between z~ and z~ is forced to be linear near the last ray through the L c term. 
Our grid has a fixed step size Az~ = h in the compactified coordinate z~, which 
corresponds to Az~ = L c /Lr h in uncompactified coordinates near z~ = 1. 

For the highest mass macroscopic black hole discussed here, M = 16, we set 
Ln = 10 9 , while for the lowest mass of M = 2~ 10 , we use Lji — 10 2 . We use 
L c = 4.096 x 10~ 9 , which can be adjusted together with Lr to obtain the desired 
resolution near the last ray. Note that Az~ w 10 _18 /i for the highest mass case; such 
a disparity in scales would have been difficult to achieve if we had used z~ as our 
coordinate even with a standard adaptive mesh refinement algorithm. We choose S to 
be between 1 and 5, the particular value of which is not too essential. 

In the + direction, for M>1, we compactify using 

z+ = Mtan(7rz+) M>1 , (14) 

with the factor of M ensuring that the singularity is not too close to the edge of the 
mesh. For M « 1, the singularity appears very close to z + = 0, so to resolve this 
region, we employ 

z+ = C z + tan p (7rz+) M < 1 , (15) 

where C z + and p are appropriate constants that again keep the singularity near the 
middle of the range of z+. For M = 2~ 10 , we use C z + = andp = 7. 

3.2. Regularization of the Fields 

It is clear from (11) that the fields diverge exponentially at 1^ and analytical re- 
sults show that they also diverge at l£ [6]. For a numerical solution then, we defined 
regularized field variables which are finite everywhere 

$ = e «( 2+ -0 (l + 4>)- M{e KZ+ -1) 

= e K( - z+ ~ z "> (l + ^ + 0o) 
6 = e K ( 2+ -0 (1 + 0) , (16) 

with (f>a = —M e KZ (1 — e~ KZ+ ). Aside from removing the divergent component 
e K ^ _z \ this definition also removes the exact classical solution M{e KZ — 1) from 
The reason for doing this came from preliminary studies that showed deviations in 
$ were small compared to the classical metric for macroscopic black holes. In terms 
of the new variables, equations (3) read 

(i + d) 2 (i + 4> + M 2 

x [d+d-4> - nd + 4> + nd-cj) - k 2 ^ + k 2 0] - Q{<j>, 6) = 

(17) 
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and 

{l + 4> + M 3 [(l + 0)d + d-6-d + 6d-6\ +Q{4>,6) = 



(18) 



with 



{ 



(i + 8) 2 [(i + 4> + fo)d + d-$ + fo)] 

(l + 9) 2 [d + (4> + 0o)d-(^ + 4> o )] (19) 

(l + 4> + 4> ) 2 [(i + e)d+d-9 - d+ed-t 



3.3. Discretization and Algebraic Manipulation 

We discretize the compactified coordinate domain as depicted in Fig. 3. A field 
a(z+, z~) is represented by a discrete mesh of values a>ij, where the indices i,j are 
integers, and related to the null coordinates through 



z c = ih < i < n p 

-A - 



= jh < j < ^ , (20) 

where h = n~ 1 is the step size in both of the compactified null coordinates. In order to 
solve the evolution equations numerically, we convert the differential equations to dif- 
ference equations by using standard, second order accurate (0(h 2 )), centered stencils: 
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(21) 

where we have introduced the notation 

f>> - d dz± d dz± f) m\ 

Once discretized, (17) and (18) give two polynomial equations which can be nu- 
merically solved for 8ij and 4>i t j, if the field values are known at the grid points 
(i, j — 1), (i — — 1, j — 1). This way, knowing the boundary conditions at 

z+ = (j = 0) an d z ~ = ~°° {i = 0)' we can calculate the field values at all points 
of the grid one by one, starting at (1,1). 



8 




Figure 3: The grid structure for the numerical calculation. We use a fixed-step-size mesh based on the 
compactified coordinates z^r , where the step sizes in both directions are equal. The emphasis on the regions 
where the fields rapidly change is attained using the compactification of the coordinates (see Fig. 2). The 
flat region before the matter pulse and the region beyond the last ray are not covered by the mesh. 



9 



Instead of solving for the two variables simultaneously (e.g. using a two dimen- 
sional Newton's method), we sum the equations (17, 18), which allows us to explicitly 
express 0jj in terms of a rational function of 6ij. We then insert this expression for 
4>i.j into (17) 1 . This way, we obtain a single variable, 10 th order polynomial equation 
for Oij. We solve this equation numerically using Newton's method, and then calcu- 
late 4>ij directly using the aforementioned rational function. Many other techniques 
are available for finding the roots of polynomials in one variable. For instance, we also 
implemented Laguerre's method, which gave similar results in terms of robustness and 
computation time. 

3.4. Richardson extrapolation with intermittent error removal 

For any function a numerically calculated on a null mesh of step size h in both 
directions, and with central differences as in (21), we have a Richardson expansion 

a h = a + c 2 h 2 + Ci h 4 + c 6 h 6 + <D(h 8 ) (23) 

where a is the exact solution, ah is the numerically obtained solution and Cj are er- 
ror functions, a, ah, and c, are all functions of z ± (we omit the explicit dependence 
for clarity), and a, Ci are independent of h. Note that we cannot prove such an ex- 
pansion exists for the class of non-linear equations we are solving, in particular if no 
assumptions on the smoothness of the initial data are made. Furthermore, we know 
the solutions generically develop singularities, thus the above series can only have a 
limited radius of convergence for generic initial data. Nevertheless, we will assume 
the expansion exists, and then, via convergence tests, check whether the solutions we 
obtain are consistent with the expansion. 

The use of second order finite difference stencils is responsible for the leading 
order quadratic convergence of the above expansion. However, using numerical solu- 
tions obtained on meshes with different discretization scales, one can obtain higher 
order convergence by using the well known Richardson extrapolation. For exam- 
ple, a fourth order convergent solution a h ,h/2 can be obtained from the following 
superposition of two approximate second order convergent solutions ah/2 an d a h ■ 
a h.h/2 = (4a/i/2 ^ o^)/3 = a + 0(/i 4 ). In theory (for sufficiently smooth solutions), 
2n-th order convergence can be obtained by an appropriate superposition of n second 
order accurate solutions, each obtained with a different mesh spacing. As we describe 
in more detail below, we use four successively finer meshes to obtain solutions that 
converge to 0(h s ) on the points of the coarsest mesh. 

Fields in the CGHS model present singular behavior, and since the position where 
the singularity first appears is a (convergent) function of the mesh size, the method of 
superposing solutions of different meshes breaks down at the first time the singularity 
appears on any of the superposed meshes. Typically, the singularity first appears on the 
coarsest mesh, and thus our domain of integration is fundamentally restricted by our 
proximity to the singularity on the coarsest mesh. Many of the physical phenomena 
we are interested in occur in this region, thus a direct use of Richardson extrapolation 



alternatively, any other independent linear combination of the equations can be used 
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for solutions over the entire computational domain does not significantly improve our 
results. To circumvent this problem, as described in more detail in the next few para- 
graphs, we instead break up the computational domain into a series of short strips in 
z~. In each strip we evolve 4 meshes, apply Richardson extrapolation to the solution 
obtained at the end of the evolution, then use this corrected solution as initial data for 
all four meshes on the next, adjacent strip. In this way the mismatch in the location of 
the singularity amongst the four resolutions is confined to be less than the size of the 
strip, which we can adjust as needed. 

Our Richardson extrapolation algorithm proceeds as follows. We divided the entire 
grid into L equal regions along z~ such that grid points i along the corresponding 
direction with j-n p < i < ^ rc p , < I < L comprise the I th region. Note that 
regions coincide at the boundaries, and here indices i and the total number of points n p 
are relative to the coarsest mesh — for finer meshes these numbers should be scaled as 
appropriate so that the I th strip spans the same coordinate volume for each resolution. 
In the I th region 

1 . We evolve the fields independently on four successively finer meshes of step size 
h,h/2, h/4 and h/8, and stop the evolution at the end of the region (i = -^n p ). 

2. At points coincident with the coarsest resolution, we calculate the appropriate 
superposition of the four meshes to give 0(h 8 ) accurate values of the fields ((f> 
and 8), and store these values on the coarsest mesh as our result. 

3. On the last (i = -^n p ) line of the region I, we also calculate the functions 
Cfe(^ ± ) to accuracy (D(h 8 ~ k ) on the coarsest mesh. We then interpolate the 
functions to the three finer meshes using successive degree four Lagrange 
interpolating polynomials. Using these interpolated values, we correct the 
field values on the finer meshes using (23). A Lagrange polynomial of degree 
4 introduces an error of order 0(h 5 ), and hence through the a term an error of 
0(h 7 ) will be introduced into the finer mesh solutions. A higher order interpo- 
lating polynomial could reduce the error, though we found that a global 0(h 7 ) 
scheme is sufficient for our purposes. 

4. We use the 0(h 7 ) accurate field values on the last (i = l -j^-n p ) line as the initial 
data for the next ((I + l) th ) region, and repeat the procedure for this region 
starting from Step 1. 

By updating the fields to more accurate values at the end of each region, the accuracy 
of the position of the singularity in the coarsest mesh is improved significantly, and the 
problem of the breakdown of the superposition near the singularity is overcome. 

We are not aware of any studies on the theoretical stability and accuracy of this 
modified Richardson extrapolation method, though our convergence and independent 
residual analysis, described in Sec. 4, shows that it works quite well, giving (for the 
most part) the expected order of convergence. 

Implementing this method, we are able to reduce the truncation error down to the 
level of round-off error using "modest" resources on a single, desktop style CPU. More 
precisely, we used 80-bit long double precision, which theoretically has a round-off 
error at the level of ~ 10~ 19 , however our Newton iteration only converged if we set 
the accuracy of the iteration to ~ 10~ 16 , which was the ultimate source of the error in 
the calculated regularized field values. When we say "round-off error" then, we will 
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mean this latter value rather than the value of ~ 10 19 one might expect from 80-bit 
precision. 

3.5. Evolution near I 

The boundary conditions on z + = are those of the vacuum and translate to 

4>{z + = 0, z~) = 6(z + = 0,z~) = 0. (24) 
For the 2^ boundary, (11) translates to 

<j>(z + ,z~ = -oo) = 9{z + , z~ = -oo) = . (25) 

The e KZ factors in the evolution equations (20) are interpreted as if e KZ is less than 
the smallest magnitude floating point number allowed by machine precision, which 
occurs for z~ < z~ rec for some z~ rec , thus, in the region z~ < z~ rec , the evolution 
equations are trivially solved by the initial conditions 6 — <f> — 0. This means it 
would make no difference if we imposed the 2^, boundary conditions on some other 
constant z~ < z~ rec line. Moreover, even if we impose the 2^, boundary conditions on 
a constant z~ line with z~ > z prec , the error introduced is exponentially small [3] and 
negligible compared to the truncation error for a certain range of \z~\. Our numerical 
method described in Sec. 3.3 sometimes fails to produce a solution for the fields in the 
early stages of the evolution near 2^ if we begin the evolution in the region z~ < z~ rec . 
We surmise the failure occurs near the line z~ = z~ rec . In such cases of failure, we 
begin the evolution at z~ ~ —5 x 10 3 , which introduces a completely negligible error. 

A related problem is that Newton's method also sometimes cannot converge to 
a solution for near 2j, even well before the last ray. Nevertheless, we were able 
to evolve the fields sufficiently close to 2^ to extract all the important asymptotic 
behavior, as described in the next sub-section. 

3.6. Asymptotic Behavior 

Many of the physical quantities of black hole evaporation are related to the asymp- 
totic behavior of $ and 6 near 2^ . The dynamical fields admit the asymptotic expan- 
sions [6] 

$ = A(z-)e KZ+ +B(z~) + 0(e- KZ+ ) 

9 = A(z-)e KZ+ +B(z-) + 0(e- KZ+ ). (26) 

which are used to calculate the affine parameter y~ on 2^ through Kexp(— ny~) = 
A(z~). The equations also admit a balance law on 2^ 

We identify the expression in square brackets to the left of the equal sign as GMb, 
with Mb identified as the Bondi mass as in [6]. We will call the term on the right hand 
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side the Ashtekar-Taveras-Varadarajan flux, Fatv ■ Asymptotic coefficients are related 
to the regularized fields through 

A(z~) = e~ KZ ~ (\ + 4>{z + = oo.O) ~ M ( 28 > 
B(z~) = lim e KZ+ (0(z+ = 00, z-)-0(z+,z- ))+M 

As mentioned in the previous section, we are not able to calculate the fields exactly 
on 1^, and evaluating the above on a line of constant z + will introduce an error of 
the order e~ KZ+ . However, it is adequate to evaluate the above at sufficiently large z + 
such that this error is less than the truncation error. It turns out the Newton iteration 
only breaks down well into the region where the truncation error dominates, and we 
calculate A on a line of constant z + in this region. 

To calculate the limit in B numerically, we need (at least) two values of z + for 
each value of z~ . From an analytical point of view, it is most desirable to use two z + 
values as large as possible. However, B is expressed as the asymptotically diverging 
factor e KZ multiplied by an asymptotically vanishing one, and calculating this via 
finite precision numerics could introduce a large round-off error due to catastrophic 
cancellation. We thus evaluate B using two z + = const lines, one the line we used 
to calculate A, the other chosen such that kz + is large enough that the fields are in the 
asymptotic region but it is also sufficiently away from the other z + = const line that 
catastrophic cancellation is not a major issue. Particular values of z + are not important. 
To calculate Mb and the ATV flux (27), we use nine-point stencils to calculate the first 
and second derivatives with respect to z~ (applying the chain rule to obtain derivatives 
with respect to z~), which have an accuracy of 0(h s ), keeping the theoretical accuracy 
of our numerical integration scheme. 

Note that since B is sub-leading relative to A in the asymptotic expansion (28), A, 
hence y~ , can be calculated more accurately. Thus, in practice we calculate the Bondi 
mass Mb by numerically integrating the ATV flux rather than directly evaluating the 
left hand side of (27). 

4. Numerical Tests 

In this section we present a few sample solutions to the CGHS model in the mean 
field approximation, and results from an array of tests we performed to ensure we are 
solving the equations correctly. 

4.1. Sample evolutions 

We calculated numerical solutions for initial black hole masses M ranging from 
2~ 10 to 16 (11). Here, we present the results for M = 8 as the macroscopic case 
for uniformity of exposition. All cases show similar convergence behavior for the 
regularized fields, though as we approach M = 16, derived physical quantities start 
to show irregular convergence patterns due to catastrophic cancellation. The fact that 
M = 8 is sufficiently large to be categorized as "macroscopic" will be established 
elsewhere, when we discuss the physical interpretation of our results. 
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The regularized fields 8 and <fi from solutions with two values of M > 1 and 
M <C 1 are shown in Figs' 4 and 5. As discussed before, a central issue with the 
numerical calculations is to ensure that we get close to the last ray and the singularity, 
as many of the interesting phenomena occur in this region. It is analytically known 
that the singularity of the CGHS model occurs when <& = y|. Moreover, $ — ^ 
evaluated on the dynamical horizon (determined by <9 + $ = 0) can be interpreted as 
the quantum corrected area of the black hole [3]. This way, we can test our proximity 
to the singularity by checking the value of the area near the singularity — see Fig.6. For 
M ;» 1, the part of our compactification scheme which emphasizes the region near 
the last ray is crucial to reach the region where the black hole area drops to less than 
a few percent of its initial value, let alone to near the Planck mass. As explicitly seen 
in the figure, had we used a uniform mesh in uncompactified z ± coordinates, a mesh 
spacing of order ft,<10~ M would have been needed. Covering a sufficient region of the 
spacetime to reveal the asymptotics would require a net coordinate range Az ± of order 
unity, implying a mesh of order 10 M points along both directions, which is of course 
impractical to achieve on contemporary computers for larger M. This important aspect 
of the problem was not clear in earlier studies, as they usually focused on M ~ 1 
[3, 4]. As we will describe in the companion paper on the physical results, the M < 1 
solutions are drastically different from the M > 1 solutions. 

A final comment — even with compactification, eventually finite precision floating 
point arithmetic will limit how large an initial mass we can simulate; with long double 
precision (80-bit), we are restricted toM<~ 20. 

4.2. Convergence of the Fields 

We compute convergence factors by comparing solutions obtained using different 
mesh spacings. Note that we are using the Richardson extrapolation scheme described 
in Sec. 3, thus in the following when we refer to a solution computed with mesh spacing 
h, h labels the coarsest resolution mesh of the four used in the numerical integration. 

First, we define 

A h f = fh - f h /2 (29) 

where fh denotes the numerical solution of a function / obtained on a grid with mesh 
spacing h. Ahf is thus an estimate, to 0(h n ), of the truncation error in /, where n is 
the rate of convergence of the algorithm. From the Richardson expansion we then get 



n = log 2 



Ah) 



log 2 



fh - fh/2 



(30) 



where the next-to-leading order term is of 0(h) because of the order of interpolating 
polynomial we use. From the above, we define an estimated convergence factor n e via 



hh — fh 
fh - fh/2 



log 2 , / (31) 



In Figs. 7 and 8 we show plots of n e for a high and low mass case respectively. An 
"issue" we have with the convergence behavior of the CGHS equations is it seems arti- 
ficially high for coarser meshes. One reason for this may be that the central difference 
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Figure 5: $ for M = 2 10 , N = 24. Left: Base-10 logarithm of <3>- j| Right: <E> — ^| at lines of constant 
z sing ~ z ~ = 10~ 4 , Again, as in Fig. 4, this shows that <E> approaches N/12, and we are 

close to the location of the last ray. Note that the field values are generally quite different from the M = 8 
case, and the singularity appears very close to 2+ = 0, which necessitated the special compactification 
scheme explained in Sec. 3.1. 
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Figure 6: Area of the black hole (<!> — yj) vs. the uncompactified distance from the last ray in a log-log 
plot for M = 8, 16 and N = 24. Note that in terms of the uncompactified coordinates, we have to be 
within Az~ ~ 10~ 8 of the last ray in order to be truly close to the singularity for M = 8, and within 
Az~ ~ 10~ 16 for M = 16. This exponential trend is general and severely limits the upper value of M we 
can use in numerical calculations if we want to reach regions "close" to the singularity. 
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Figure 7: Convergence of the M = 8, N = 24 case: n e (z ± ) for h = (left) is mostly in the range 

9 — 10, and for h = (middle) is around 8. For h = 2~ 12 (right) we reach machine round-off, and 

thus loose convergence, hence the "noisy" pattern. 
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Figure 8: Convergence of the M = 2~ 10 , N = 24 case: n^z*) for h = 2 -10 (left) is around 10, and for 
h = 2 -11 (middle) is around 8. Again, as with the M = 8 case in Fig. 7, for h = 2~ 12 (right) machine 
round-off error begins to dominate the error, hence the "noisy" pattern. This effect is already visible in 
certain regions of the h = 2 -11 case. For lower mass black holes, round-off is reached with coarser meshes 
relative to the higher mass black holes. 
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scheme (21) we use solves the homogeneous part of the wave equation (d+d-f = 0) 
exactly (to within round-off), irrespective of the step size. Furthermore, with our 
choice of variables and regularization scheme, it is only the non-linear quantum correc- 
tions that introduce non-trivial evolution, and initially the effects of this will be small. 
Though regardless, in the limit of zero h we should approach the expected convergence 
behavior; as shown in these figures, we do see this trend, though have not quite reached 
the limiting behavior before machine round-off error is reached. 

As mentioned, reasons for the anomalous convergence behavior may be the com- 
pactification and special initial data we choose, namely regularized fields that are ini- 
tially adapted to the classical solution. To check this, we evolved a test case where we 
imposed the initial conditions for M = 11, N = 11 at z+ = 0.25 rather than X^. Note 
that this is not a physically correct solution as it will violate the constraints, though 
it is mathematically perfectly valid non-trivial initial data for the evolution equations. 
We set the domain of computation to z~ e [0.25, 0.5] and z+ G [0.0, 0.25] to avoid 
any singular behavior. Using four meshes for the Richardson extrapolation in this test, 
the truncation error was again reduced down to round-off level for even the coarser 
meshes, so for this test alone, we only employed three successively finer meshes in the 
extrapolation scheme; hence the resulting truncation error is expected to scale as h 6 . 
The result is shown in Fig. 9, where we see the expected convergence. We also tested 
two-mesh Richardson extrapolation for the same case, and obtained the expected h A 
convergence. 

4.3. Convergence of Physical Quantities on 

The physical quantities we are interested in, including y~ (z~), Fatv an d Mb, are 
all functions of the fields, thus in theory they should inherit the convergence behavior 
of the fields. Some of these quantities require computing first and second derivatives 
of the fields, and so to maintain the theoretical convergence factor of 7, one should 
use 9-point finite difference stencils. However, catastrophic cancellation plagues the 
numerical derivatives near the last ray, as the regularized fields vary extremely slowly 
in this region, and this seems to be the limiting factor in the accuracy in which we 
can compute physical quantities. Though in general we do not need high order conver- 
gence of derived quantities to achieve high accuracy. A case-and-point is Mb, obtained 
by integrating Fatv ■ Fatv is dominated by round-off near the last ray in most cases, 
though, once integrated over y~, this region contributes insignificantly to Mb- Further- 
more, simple trapezoidal integration is adequate to achieve quite accurate estimates of 
M B , as illustrated in Fig. 10. 

4.4. Independent Residuals 

As a final test of the code, from numerical solutions, we compute independent 
residuals of the differential equations (17) and (18). Specifically, we calculate the 
derivatives using three-point stencils centered at the mesh points, rather than the cell 
centered differences used for the solution. Three- point stencils limit the convergence 
of the independent residual to quadratic order, regardless of the convergence of the 
fields themselves. We observe the expected quadratic convergence in all cases. 
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Figure 9: The convergence factor n e for h = 2 8 where as a test we imposed the (unphysical) initial 
conditions for M = 11, N = 11 at = 0.25 rather than J^. We only evolved the fields in the region 
Zc 6 [0.25, 0.5], 2+ £ [0.0,0.25]. This solution is not physically relevant, though tests the behavior of 
the numerical code away from any of the null infinities or singularities. Here, for each base resolution, three 
meshes where used in the Richardson extrapolation scheme, which should give 0(h 6 ) convergence, and 
does to good approximation as shown in the figure. 
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Figure 10: | j for various values of h for M = 2 10 (left) and M = 8 (right). For most of the range, 

there is clear quadratic convergence. The dominant error here is from the trapezoidal integration method, 
and not a reflection of the truncation error from the numerical calculation of the fields. 
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5. Conclusions 

In this paper we presented some details of the technical aspects of our numerical so- 
lution of the CGHS model. Beside the links this model provides to black hole evolution 
in 3 + 1 dimensions (which will be published elsewhere), it also presents considerable 
numerical challenges to solve, in particular in the macroscopic mass limit. The fact 
that the CGHS model has two distinct regimes in parameter space, M > N/24 and 
M -C TV/24, has not been emphasized in the literature before. Not only do these two 
regimes have radically different physical properties and interpretations, their numeri- 
cal analysis also presents considerably different levels of challenge. Existing numerical 
studies of the CGHS model focused on the intermediate mass range M ~ N/24, for 
example ^fjy = 1 in [3] and = 2.5 in [4]. This case, like the low mass region, is 
considerably easier to solve numerically. In this regime, most of the evolution of the 
fields is not confined to a small region near the last ray, thus, the uncompactified z ± 
coordinates are adequate to cover the quantum-corrected spacetime. However, many 
of the interesting phenomena of the black hole evaporation cannot be observed in this 
regime. In the physically more interesting case of macroscopic black holes, we are 
led to using compactified coordinates: we need to start the calculation sufficiently far 
away from the last ray, yet at the same time have high resolution (Az~ ~ 1CT 16 for 
M = 16, N = 24) near the last ray. 

A correct estimate of the position of the singularity necessitates very low trun- 
cation error, which we obtain using Richardson extrapolation (with intermittent error 
removal). This takes the results from four successively finer meshes to obtain results 
that theoretically scale as 0(h 7 ), and we use this scheme to reduce the truncation error 
to the level of machine precision (10~ 16 ). Even though there is still space for im- 
provement of our method using improved compactification schemes, higher precision 
numerics, etc. we have enough accuracy to, for the first time, discern the physics of the 
CGHS model near the last ray in the macroscopic mass regime. 

Acknowledgments: We would like to thank Abhay Ashtekar who suggested the 
problem to us, collaborated with us on the physical aspects of the CGHS model, and 
made useful suggestions on this manuscript. We would also like to thank Amos Ori for 
useful discussions. We acknowledge support from NSF grant PHY-0745779, and the 
Alfred P. Sloan Foundation (FP). Some of the simulations were run on the Woodhen 
cluster at Princeton University. 

References 

[1] C. G. Callan, S. B. Giddings, J.A. Harvey and A. Strominger, Phys. Rev. D 45, 
R1005 (1992). 

[2] S. B. Giddings and A. Strominger, Phys. Rev. D 46, 627 (1992); Phys. Rev. D 47, 
2454(1993);J.G. Russo, L. Susskind and L. Thorlacius, "Black Hole Evaporation 
in 1 + 1 Dimensions", Phys. Rev. D 46, 3444 (1992); Phys. Lett. B 292, 13 (1992); 
L. Susskind and L. Thorlacius, Nucl. Phys. B 382, 123 (1992); A. Strominger, 
Phys. Rev. D 46, 4396 (1992); S. P. de Alwis, "Black hole physics from Liouville 
theory" Phys. Lett. B 300, 330 (1993); Phys. Rev. D 46, 5429 (1992); A. Bilal 



23 



and C. Callan, Nucl. Phys. D 394, 73 (1993); S. W. Hawking and J. M. Stewart, 
Nucl. Phys. B 400, 393 (1993); 

[3] T. Piran and A. Strominger, Phys. Rev. D 48, 4729 (1993). 

[4] D. A. Lowe, Phys. Rev. D 47, 2446 (1993); 

[5] S.B. Giddings, arXiv:Hep-th/9412138; A. Strominger, arXiv:Hep-th/9501071. 

[6] A. Ashtekar, V. Taveras and M. Varadarajan, Phys. Rev. Lett. 100, 21 1302(2008) 

[7] J. Winicour, Living Rev. Rel. 1, 5 (1998) [arXiv:gr-qc/0102085]. 



24 



